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Abstract 

A numerical method for computation of heteroepitaxial growth in the presence of 
strain is presented. The model used is based on a solid-on-solid model with a cubic 
lattice. Elastic effects are incorporated using a ball and spring type model. The 
growing film is evolved using Kinetic Monte Carlo (KMC) and it is assumed that the 
film is in mechanical equilibrium. The strain field in the substrate is computed by 
an exact solution which is efficiently evaluated using the fast Fourier transform. The 
strain field in the growing film is computed directly. The resulting coupled system is 
solved iteratively using the conjugate gradient method. Finally we introduce various 
approximations in the implementation of KMC to improve the computation speed. 
Numerical results show that layer-by-layer growth is unstable if the misfit is large 
enough resulting in the formation of three dimensional islands. 

1 Introduction 

Epitaxial growth is the process where crystals are grown by the deposition of atoms 
in a vacuum. Typically the deposition rate is small and the crystal is grown, loosely 
speaking, one layer at a time. In this paper, we consider the computation of strained 
epitaxial growth when the strain arises because the natural lattice spacing of the sub- 
strate and the deposited material are different. This difference is called the mismatch. 
Heteroepitaxial growth is experimentally observed to grow in the following growth 
modes: 

1. Frank- Van der Meer growth: crystal surface remains fairly flat, growth occurs in 
the layer-by-layer fashion. 

2. Volmer- Weber growth: three dimensional islands form on the substrate without 
a wetting layer. 

3. Stranski-Krastanov growth: the film grows in a layer-by- layer fashion for a few 
layers, and then Volmer- Weber growth begins. This results in three dimensional 
islands on top of a wetting layer. 
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Figure 1: Germanium on Silicon- Due to elastic interaction. The bottom configuration has 
less energy than the top one 

The type of growth mode depends on many parameters, an important one is the 
mismatch. In many cases, when the mismatch is high, one finds Volmer- Weber growth 
and when the mismatch is small layer-by-layer growth is observed. For intermediate 
values of the mismatch Stranski-Krastanov growth is often seen. For an overview see, 
for example, Ref. PH ITlj. 

In homoepitaxy, the effects of strain are usually very small and quite often ignored 
in the many models. In general, the morphology of a growing film by homoepitaxy is 
reasonably well understood. It is known that in some cases a homoepitaxially grown 
film can undergo an instability resulting in mound formation. Typically these phenom- 
ena are due to kinetic effects, for example a step-edge barrier (21 El El or enhanced edge 
diffusion j£] can cause mound formation. 

However, when a species of atoms grows on a substrate of a different species, for- 
mation of 3D islands is observed in many situations. It is generally believed in many 
cases (for example for the growth of Germanium on Silicon) that this is a thermody- 
namical effect. In particular, the elastic energy stored in a strained flat interface is 
greater than when there are three dimensional islands. This is due to the fact that 
in the latter case the atoms have more opportunity to relax (see Figure However, 
the surface energy of three dimensional islands is greater than that of a flat interface. 
This implies that the morphology of heteroepitaxially grown films is determined by the 
interplay between elastic energy which is a bulk effect and surface energy which arises 
from broken bonds. 

1.1 Modeling Elastic Effects 

Elastic effects in thin films can be studied with fully continuum models or Burton- 
Cabrera- Frank [2] type models that consider elastic effects between steps. In this paper, 
we shall consider a fully discrete model which is evolved in time using a kinetic Monte 
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Carlo method. Naturally, such an approach has the disadvantage of not being able to 
simulate on large length scales. However, it offers the advantage that nanoscale physical 
effects such as island shape fluctuations and nucleation are naturally incorporated. 
One of the first, if not the first, model in this direction is due to Orr et al. JB]. They 
accounted for the elastic interactions using a ball and spring model, which takes into 
account nearest neighbor and next nearest neighbor interactions. This was combined 
with a solid-on-solid type model which was then used with KMC to simulate a growing 
heteroepitaxial film in 1+1 dimensions. If the misfit was below a critical value, the 
film grew in a layer-by-layer fashion. On the other hand, if the misfit was above the 
same critical value, then the film was observed to grow in the Volmer- Weber fashion. 
Later Lam, Lee and Sander [Hj provided a more efficient implementation of this model, 
which allowed them to perform simulations using parameter values that were more 
physically reasonable and compute for larger domains. This work has been recently 
been extended to three dimensions [S]. 

Ratsch et al. ^Hj studied three dimensional heteroepitaxy, however they did not 
take explicitly into account the harmonic forces between atoms, but rather they used 
an approximate treatment ^7] based on the Frenkel-Kontorova model. The model was 
used to investigate the island size distribution in heteroepitaxial growth [TH] . 

Off lattice KMC simulations of heteroepitaxial growth in 1+1 dimensions where pre- 
sented in a series of papers ^3 El El- In these computations the forces between atoms 
were modelled using Lennard-Jones interactions. The misfit is easily incorporated by 
changing parameters in the potential. One advantage of this approach is that dislo- 
cations are naturally included, which is not the case with the ball and spring model. 
These simulations also demonstrate that if the misfit is sufficiently large, layer-by-layer 
growth is unstable and mounds form. 

A more sophisticated discrete elastic model was introduced by Schindler et al. j2Qj. 
This model is based on a discrete form of the continuum elasticity equations. The 
approach presented here could be used to solve their model as well. 

2 Model description and Kinetic Monte Carlo 

The model we shall use is a three dimensional version of the model proposed in Refs. 
|13l For the convenience of the reader we shall now describe this model. To fix 
ideas we shall assume that the deposited atoms are Germanium and the substrate is 
composed of Silicon. The atoms occupy sites arranged on a simple cubic lattice with 
no over hanging atoms allowed. This means that the height of the surface is a function 
of the substrate location. We assume that atoms bond with their nearest and next 
nearest neighbors. Each atom can be linked to its six nearest neighbors located at a 
distance a, and to its twelve next neighbors located at a distance a\/2. For example, 
an atom on a flat plane surface orthogonal to one of the coordinate axis will have five 
bonds with nearest neighbors, and eight bonds with next nearest neighbors, while an 
atom sitting on top of that same flat surface will have five bonds (one with a nearest 
neighbor and four with next near neighbors). We shall assume the chemical energy 
associated to all these bonds is the same. The total chemical bond energy associated 
to each atom is therefore Ef, = —^Nj,, where Nb is the number of bonds of each atom, 
and 7 the energy associated to each bond. 

The elastic effects in this model are taken into account by assuming that the bonds 
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will act like a spring between the atoms. We will use a s and a g to denote the lat- 
tice spacing between Silicon and Germanium atoms respectively. We shall denote 
respectively by ki, and ko the spring constants corresponding to longitudinal (near- 
est neighbor) and diagonal (next nearest neighbor) bonds. For ease of exposition, we 
shall assume that both Silicon and Germanium have the same spring constants. Since 
o-g 7^ a>s mechanical, force will arises (the calculation of which is described in detail be- 
low). In many systems, the time taken for sound waves to propagate across the sample 
is much smaller than the time scales associated to any growth process. Therefore we 
assume that our mass-spring system is always in mechanical equilibrium. 
Each atom, p, of the system will hop with a rate T m [Sj given by 

r m = #0 exp (l^r) (!) 

where 

AE = ^(without atom p) — .E(with atom p) (2) 

is the change in energy of the entire system when atom p is completely removed. Ro is 
the attempt frequency, ks is the Boltzmann factor, and T is the lattice temperature. 
Since the chemical bonds are purely local, then we can write (J2J) as 

AE = n b 7 - AE e i as 

where rib is the number of chemical bonds of the atom, 7 is the energy associated to 
the chemical bond, and 

AE e i as = E e i as (with atom p) - E eias (without atom p) (3) 

We note that AE e i as is always nonnegative and when combined with © implies that 
elastic effects will always increase the hopping rate. 

We shall evolve the model in time by the use of kinetic Monte Carlo (KMC). The 
basic KMC method can be described as follows. 

1. Pick an atom at random on the surface, uniformly among all surface atoms 

2. Compute its number of bonds 

3. Compute the contribution of the elastic energy associated to the atom 

4. Pick a random number, r, uniformly distributed in [0, 1]. 

5. If r < exp(—(ribj—AE e i as —K)/kBT) then move the atom, according to a random 
direction, chosen uniformly among all possible directions 

The constant K is computed in such a way that the numerator rib'j — AE e [ as — K 
that appears in the exponential is always non negative. This sets the fastest rate in 
the problem to unity. The method described is based on rejection. It is well known 
that if the number of possible event types is small (as is the case for this model when 
there are no elastic effects) then rejection- free Kinetic Monte Carlo can provide a much 
more efficient algorithm. 

While this model is idealized, it nevertheless captures the essential physicals effects 
of heteroepitaxial growth, such as adatom diffusion, nucleation, surface diffusion, long 
range elastic interaction. In addition, since the model is evolved in time using kinetic 
Monte Carlo, it naturally captures effects associated with fluctuations. 
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Figure 2: The reference configuration is obtained by compressing the Germanium atoms to 
have the same horizontal spacing as the Silicon atoms. The vertical spacing is chosen so that 
Germanium atoms are in equilibrium. 



3 Elastic computations 

The main difficulty in the implementation of this model is the computation of the 
strain field. In this section we shall outline our approach for solving this problem. For 
the basic set up we follow Lam et al. |S], as is described in next section. However, our 
numerical implementation is different from the one used in jS]. One important feature 
of our work is that we provide an exact solution for the elastic displacement in the 
substrate which is efficiently evaluated using fast Fourier transforms. 



3.1 The reference configuration 

The reference configuration we choose consists of a periodic array of complete layers 
of Germanium atoms on top of a periodic array of Silicon. The Germanium atoms are 
compressed so that their horizontal lattice spacing matches that of the Silicon atoms, 
see Figure |31 The vertical lattice spacing, , is chosen so that the resulting system 
is in mechanical equilibrium. We will now describe the computation of ai in two 
dimensions. It is useful to introduce the following dimensionless quantity 

a g - a s 
e = — . 

a g 

which is denoted as the mismatch. Typical values of e range from -0.05 to 0.05. For 
example the mismatch for Germanium on a Silicon substrate is 0.04. In order to 
deduce the atom displacement with respect to the reference configuration we need to 
compute the forces experienced by an atom due to each of its neighbors. Elementary 
considerations allow to compute that, to first order in the ratio e, one has 

F 1 = F h (qJ,F 2 = F dv (^),F 3 = F v (® 1 \F i = F DV (~^\f 5 = -F lt 

where F v = k L (a L - a s ), F H = k L (a g - a s ), and F DV = k D (2a g - a L - a s )/2. 

The value of ai is determined by requiring that these five forces sum to zero for 
atoms in the reference configuration. By symmetry, the forces in the x direction sum to 
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Figure 3: The reference configuration is obtained by compressing the Germanium atoms to 
have the same horizontal spacing as the Silicon atoms. The vertical spacing is chosen so that 
complete layers of Germanium are in equilibrium, b is the net force on an atom due to the 
compression. In principle b can be nonzero for any Germanium or Silicon atoms in the top 
row, if the top layer is not complete 



zero. On the other hand balancing the z components of the force one has 2Fjjy = Fy 
which implies 

k D (2a g - a L - a s ) + k L (a g - a L ) = 0, 



which gives the following expression for ql 



aL 



1 + e 



k D 



k L + k D 



A similar argument can be applied to the three dimensional lattice. In this case 
each atom can interact with six nearest neighbors, located at a distance a, and twelve 
diagonal next to nearest neighbors, located at distance ay/2. The interaction with the 
8 corner neighbors, located at a distance ay/3, are neglected. 

As in the two dimensional case, we shall denote by ki and ko the two spring 
constants corresponding to the interaction between nearest neighbors and diagonal 
neighbors. Each bulk atom is surrounded by 18 neighbors, (six longitudinal and twelve 
diagonal). We denote by the contribution of the force on a given atom due to the 
presence of its neighbor in the direction (i,j,k). For example, the 3D equivalent of 
force F3 of Figure 01 would be Fo,o,-i- 

The six forces aligned along the coordinate axis have the expression 



ijk 




with k G {—1, 0, 1}, where \i\ + |j| + \k\ = 1. (4) 
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Fy and Fh are given above. The twelve diagonal forces are given by 



and 




FiOk = -\ , with t,fc€{-l,l}, (5) 



bjfc = - jFdv , with i, fee {-1,1}, (6) 



Fijo = -\ jF DH , with y-G{-l,l}, (7) 



Fdv is given above and Fdh = ko{a g — a s ). It is convenient to set F^ = if 
|i| + |fc| = 3. 

As in the two dimensional case these forces must sum to zero in the reference 
configuration. The x and y components will vanish by symmetry. The forces in the z 
direction vanish if 

F v + 4F DV = 0, 

which implies 

2k D (2a g — a,L — a s ) + ki{a g - a£) = 0. 
This gives the following expression for cll, 

ol = "»( 1+t ^fk)- (8) 

3.2 Computation of the interaction 

Let us denote the displacement, with respect to the reference configuration, of an atom 
at site k) by the vector (uejk, V£jk,W£jk) and the force experienced by this atom as 
{fijk-, 9ijki htjk)- This force will arise from the interaction of the atom with its nearest 
neighbors and next nearest neighbors. For example the x component of the force is 
given by 



fijk = k L ([ue+ijk - u£ jk ] + [ug-ijk - u£j k ]) 
ku 

+ ~y {[ue+ijk+i - utjk] + [ue-ijk+i - u£jk\) 

+ {[ue+ijk-i - uejk] + [ue-ijk-i - uijk}) 
ko 

+ ~y {[ue+ij+ik - uejk] + [ue-ij+u - uejk]) 
ko 

+ ~y~ (N+ij-ifc - uejk] + [ue-ij-ik - u£jk]) 

+ (N+lj+lfe ~ v ijk] + [Vi-lj-lk - V£jk\) 

~ ~Tr (N+lj-lfe _ v ijk\ + [ve-lj+lk - v ijk\) 



+ ~2 ([ w e+ijk+i ~ wejk] + [w£-ijk-i - wijk]) 
- -y ([we+ijk-i ~ W£ jk ] + [we-ijk+i - U£ jk \) 

~\~ ^ Fmnq ' &x (9) 

(m,n,q) £neigh(£ ,j, k) 

where F mnq are given by Q4l7j) . Each term in square brackets and each F mnq represents 
the interaction of an atom at site (£, j, k) with potential nearest and next nearest 
neighbors. If no such neighbor exists then the term should not be included. 

Suppose we have N atoms and we denote the relative displacement of the pth atom 
by Up and let F p denote the force it experiences. We also let b p denote the sum of all 
forces given by (@EJ), acting on the atom when its position is the reference configuration. 
Next we define the following vectors in M 3Ar , u = (u\, . . . , un) t , b and F are similarly 
defined. Then we can write 

F = ^u + b. 

We remark that for atoms that are completely surrounded by other 18 atoms, or atoms 
that are on a horizontal surface, the corresponding b is zero, since all the forces acting 
on them sum up to zero, which is consistent to the fact that a rectangular box of atoms 
in the reference configuration is in equilibrium. As a consequence, the vector b has 
nonzero elements only for atoms at the surface. The matrix vector product, ^4u can 
be deduced from © and similar relations for gijk and htjk- 

The equilibrium position of atoms in a given configuration is obtained by setting 
F = 0, i.e. by solving the large linear system 

„4u + b = 0. 



3.3 Contribution of the substrate 

It is known that elastic interactions can be very long ranged. For example the elastic 
interaction between two island behaves like d~ 2 where d is the distance between the 
island centers ^H]- This indicates that elastic interaction can penetrate deep into the 
substrate. On the other hand, the interaction range is certainly much shorter than 
the thickness of the substrate. For this reason it is prudent to consider the substrate 
to be semi-infinite in the z-direction. To reduce boundary effects we consider periodic 
boundary conditions in both the x and y directions. In this section we shall derive a 
formula that expresses the force on the surface atoms of the substrate completely in 
terms of their displacement. 

The surface of the substrate corresponds to k = 0, and the atoms of the bulk 
substrate will be indexed using negative k values. In the substrate (k < — 1) all atoms 
have a complete set of neighbors, consequently we can write the force, in component 
form, on the atom at site (£,j, k), k < —1, as 

fijk = ki(ue + ijk — 2u£jk + ue-ijk) 
ko , 

"I — 2~ Wfijfc+i + u £-ijk+i + u e+ijk-i + ui-ijk-i 
+ u e+ij+ik + ut-ij+ik + u e+ij-ik + ue-ij-ik - Suejk) 
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+-^-( v e+ij+ik + ve-ij-ik - ve+ij-ik - ^-ij+ifc) 

+-^-{we+ijk+i + we-ijk-i - we+ijk-i - we-ijk+i), (10) 



gtjk = k L (v£j + i k - 2v£jk + Vij-ik) 

"I — 2~( w #+ifc+i ~*~ + v £j+ik-i + 

+^£+ij+ifc + V£-ij+ik + vi+ij-ik + vi-ij-ik - 8vejk) 

+-^-( u e+ij+ik + ue-ij-ik - ue+ij-ik - w-ij+ik) 
+-^-(wej+ik+i + wtj-ik-i - wej+ik-i - wej-ik+i), (11) 

fyjfc = k L (wejk+i - 2wejk + wgjk-i) 

+yWiHi + ^i+ifc-i + w tj-\k+i + wej-ik-i 
+W£ + ij k+1 + wg-ijk+i + w^+ijfc-i + wg-ijk-i - Swijk) 

^ — ^-( u e+ijk+i + ug-ijk-i — u £+ijk-i - ug-ijk+i) 

+ ~^-(vej+lk+l + V£j-lk-l - V£j+lk-l - V£j-lk+l)- (12) 

At the surface of the substrate (k = 0), one has a slightly different expression, because 
there are no atoms on top, namely: 

fejo = kL(u£+i t jfi — 2 u e,j,o + ue-i,j,o) 
kD i 

{ u e+i,j,-i + u£-i,j,-i + ue+i,j+i,o + ue-i,j+i,o 
+u£+i,j-i,o + ue-ij-ifl - 6u£ :j:0 ) 

+ ^-( V i+l,j + l,0 + V£-lj- lfi - V£ +1>j -l >0 - V£-lj +1>0 ) 

+ ^{W£-I,j,-1 - W l+1 j-x), (13) 

9£jo = k L (v£ ij+lfi - 2v e j fi + V£j-i,o) 
kD , 

+~2" W+ij+i,o + ve-ij+1,0 + V£+ij-ifl + V£-ij-ifl 
+ve,j+i-i + V£j-i-i - 6v£ jjt o) 

+-^-( u e+i,j+i,o + U£-ij-ifi - U£ + ij-i t o - U£-ij+i j0 ) 
kD 

+yhj-i,-i " (I 4 ) 

hjo = k L (w£ tj _i - W£ Jfi ) 

+-^{we+i,j,-i + m-ij-i + u^j+i _i + w£j-i-i - 4:we jjj0 ) 



+-^{ut-i,j,-i - ue+i,j,-i) 

+-7r( v e,j-i,-i - V£,j+i,-i)- (15) 

Let us now consider a Fourier expansion of the displacement in the x and y direction. 
The generic Fourier mode will take the form 

u tjk = u k (^v)e m+jr >\ 

Vijk = Vk&tfeW+M, (16) 
wtjk = Mt,v)e i( # +i,,) . 

By inserting this Fourier expansion in the expression of the surface force (|13I15|I ° ne 
obtains the relations 

fo = 2/c^uo(cos £ — 1) + ko[u-i cos £ + uq{2 cos £ cos 77 — 3) 
— 2vq sin£ sin 77 — iw-\ sin£], 

go = 2kLVo(cosrj — 1) + ku\v-\ cos 77 + wo (2 cos £ cos 77 — 3) 

— 2u sin £ sin 77 — iw_i sin 77], (17) 

ho = kiiw-i - wo) + k D [w-i(cos^ + cos 7/) - 2w ] 
— i(u- \ sin£ +V-i sin 77)]. 

In the relation above we have suppressed the dependence of all Fourier modes on (£, 77). 

Ea.(|17|) gives a relation between the Fourier modes of the force and the Fourier 
modes of the displacement. However, our goal is to express (/o,<?o>^o) m terms of 
(uo,vo,wo). Once this is done then the force field at the surface can be computed 
from its Fourier modes by inverse discrete Fourier transform. In order to accomplish 
our goal, we need to express u-x,v—i,W—i in terms of uq, vq, u>o, and substitute their 
expression into (|17jl. This can be done as follows. First, let us insert the Fourier 
expansion ()16l) into ()1UH12|) obtaining: 



f k = 2k L u k (cosi-l) + k D [(u k+ i + u k _ l )cosi + u k {2cosr]cosi-'i)] 
+ikf)(w k+ i — uik-i) sin£ — 2kf)V k sin£ sin 77 

g k = 2k L v k (cosr] - 1) + k D [(v k+1 + v k _i) cost? + % (2 cos 77 cos £ - 4)] (18) 
+ikD{w k+ \ — u>fc-i) sin 77 — 2kr>u k sin£ sin 77 

h k = k L (w k+ i - 2w k + w k -i) 

+kD[(w k+ i + 7i>fc-i)(cos£ + cos 77) - 4w k ] 
+ik D [(u k+1 - Ujt_i) sin£ + (v k+1 - v k -i)smr]] 

The discrete equations given by (|T8|) are solved using the following substitution 

u k = ua k , v k = va k , w k = wa k , (19) 
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where we look for solutions with \a\ > 1, since we expect the Fourier modes to decay 
as k —* — oo. Inserting this ansatz into the expression ([18)1 of the Fourier modes of the 
force acting on the inner points of the substrate, one obtains: 



(20) 











| = n(a) | 











where the entries of the matrix $7 are given by 

u>n = 2&l(cos£ — 1)« + /c£)[cos^(l + a 2 ) + 2(cost7cos£ — 2)a] 

0J22 = 2/cl(cos rj — l)a + ko[cos rj(l + a 2 ) + 2(cos rj cos £ — 2)a] 

^33 = M" 2 - 2a + 1) + k D [(a 2 + l)(cos£ + cos??) - 4] 

Ui2 = t°2i = — 2a/co sin ^ sin r/ 

^13 = ^31 = ik D (a 2 - 1) sin£ 

^23 = ^32 = ikoia 2 - 1) sin r/ 



Note that matrix 0, is symmetric, but not self-adjoint. 

Since all forces in the bulk have to be zero (all such atoms are in mechanical 
equilibrium), then one has 

/n\ 

n\ v =o. (21) 

This homogeneous system has nontrivial solutions only if 

P(a) = det(fi) = 0. (22) 

This relation results in an algebraic equation for the values of a. The polynomial P{a) 
is a six degree polynomial, therefore it admits, in general, six roots. Note that, because 
of the structure of the matrix Q, matrix a 2 Q.{l/a) is equal to Q(a) with W13 = Usi and 
W23 = W32 of opposite sign. This does not change the expression of the determinant, 
and therefore if a ^ is a root, then also 1/a is a root. This means that the number 
of roots a such that \a\ > 1 is equal to the number of (nonzero) roots a such that 
\a\ < 1. The roots that are of interest for us are the ones that decay as k — > oo, i.e. a 
: \a\ > 1. 



3.3.1 Eigenvector computation 

In this subsection we describe a general procedure for the computation of the eigenval- 
ues and eigenvectors, that works also in the case of multiple eigenvalues. 
The goal is to solve the problem given by (|21j) . which we write as 

S> = (23) 

where r is a three-component vector (we drop the arrow on top), and to find inde- 
pendent eigenvectors even if some eigenvalues coincide. First compute the eigenvalues 
by solving the algebraic equation (|2*2*|) . Consider the three eigenvalues ai such that 
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\ote\ > 1, £ = 1, ... ,3. If they are all distinct, then the three eigenvectors correspond- 
ing to them will be independent. If two of them are coincident, let us say 02 = 03, 
then one has to find two independent eigenvectors corresponding to the coincident 
eigenvalues. 

A unified treatment of the problem is obtained by the use of the singular value 
decomposition (SVD) of matrix f2. The procedure works as follows. First compute 
a±, 02,03. If they are distinct, for each of them compute &e = fi(a^), £ = 1, ... ,3. 
Perform the SVD of Qf. = WEV', where U and V are unitary matrices (i.e. 

is a diagonal matrix containing the singular values of fig. 
Taking into account that U is non singular, problem Q23|) reads 

evV = 0. 

Since is singular, then £ = diag(oi, 02, 0), therefore one has 

ai v tr )i = °' ff2 v tr ) 2 = °' (^ tr ) 3 = arbitrar y 

Let us choose (V 'r) 3 = 1. 

Assuming a<i 7^ 0, i.e. that the matrix has rank 2, then one has 



therefore 




i.e. r is the third column of V. 

If two roots are coincident, say 01,02 = 03, then first compute the eigenvector r\ 
using the procedure above applied to matrix fi(ai). For the computation of the other 
eigenvectors there are two possibilities: either the rank of the matrix = ^3 is 1, 
i.e. o"2 = 0, or the rank of the matrix is 2, i.e. 02 7^ 0. However, the latter case 
never happened in all our computations, and we conjecture it can never happen for 
our problem. Therefore we assume that 02 = 0. Repeating the procedure above, one 
finds that T2 and r$ can be computed, respectively, as the second and third column of 
the matrix V . 

3.3.2 Surface Force Formula 

We denote by a solution of the system 

n(a e y=o, 

with £ = 1,2,3. Then we can decompose the vector (uo,vo,wo) T on the basis of the 
eigenvectors, i.e. 

^0 \ 

VQ = Clfi + C 2 f 2 + C 3 T3. 
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Once the constants c\, 02,03 are computed, one can write 




Cl _ . C 2 _ C 3 _ 

= — n H r 2 H r 3 . 

ai 02 «3 

The two relations allow to express u~i, v~i, w~i in terms of uo, vo, wo- Once this is 
done, one can substitute this expression into (|17[) and obtain the final relation between 
(/o) 30)^0 ) an d (uo,vo,u>o)- This relation has to be computed for all Fourier modes 
(£, 77). Periodicity implies that £ = 2-7rm/M, = 2irn/M, m, n = 1, . . . , M. 

The complete algorithm for the computation of /«o from u^o can be summarized 
as follows. 

Computation of fgjo from it^o 

0. Preprocessing. Given M, for each mode (1711,1712), Tn\, m 2 = 1, . . . , M, solve the 
eigenvalue problem (|21j) . and store the eigenvalues and eigenvectors. 

1. Given u^o, perform the discrete Fourier transform in I and j and compute all 
Fourier modes uo,vo,wo. 

2. For each mode, compute U-i,v-\, using pre-computed values of eigenvalues 
and eigenvectors. 

3. Compute the Fourier modes of the force f,g,h, using Eq.(|17|) 

4. Compute the force by inverse discrete Fourier transform 

All discrete Fourier transforms can be efficiently computed by FFT algorithms in 
0(M 2 log M) operations. In all our calculations we used the FFTW package devel- 
oped at MIT p. 

3.4 Elastic Displacement Computation 

Let us assume that we have deposited iV atoms on a substrate of size M x M. Let 
us use u g £ M? N to denote the relative displacement of the Germanium atoms on the 
substrate. We use u s S M 3A:f2 to denote the relative displacement of the top layer of 
atoms on the substrate. Then the equilibrium position of the particles can be obtained 
by solving the following linear system 

:;)-(^2)(iH5:)- 

The matrices appearing in the system have the following meaning. The forces acting 
on the M 2 Silicon atoms on the surface of the substrate are 

F s = Su s + Bu g + b s . 

Here Su s is the force on the atoms at the surface of the substrate due to all the (Silicon) 
atoms in the substrate. This is efficiently computed using the results from the previous 
section. Bu g is the force on the substrate surface due to the Germanium atoms, and 
b s is the sum of the forces given by (|4I7|) . The force acting on the N Germanium 
atoms on the substrate is given by 

Fg = B T U S + AUg + bg, 
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where Au g are the forces that arise from the interactions between the Germanium 
atoms, B T u s is the force on the Germanium atoms due to the top layer of Silicon 
atoms, and b 9 is the sum of the forces given by (|4I7|) . 
We observe that the matrix 

S B 



B T A 



(25) 



is a symmetric negative semi-definite matrix; it has 3 zero eigenvalues, corresponding 
to the free translation in the 3 directions of the coordinate axis. The system is clearly 
invariant for translation along the directions parallel to the substrate. It is also invari- 
ant along the direction orthogonal to the substrate, because the substrate is considered 
semi-infinite. This can be understood by the following argument. For a substrate of a 
finite thickness, let us say of Nl layers, a unit displacement in the direction orthogonal 
to the substrate will produce an elastic force per unit atom equal to 

F = ^-{k L + 2k D ) 

which vanishes as Nl — > oo. Therefore no resistance is opposed to any translation. 

Notice that the matrix B T is the transpose of matrix B G M 3M x3iV . A is a 3N x 3iV 
matrix. A and B are sparse and the matrix vector products are efficiently evaluated 
using expressions similar to ©• System (|24j) can be solved by an iterative scheme 
for large, sparse linear systems, making use of the symmetry and definiteness of the 
coefficient matrix. Here we shall use just a simple conjugate gradient method, leaving 
the search for a more efficient method to future investigations. 



4 Evaluation of the elastic energy 

Once the strain field is determined, the elastic energy is computed as follows. The 
energy associated to the bonds is given by 

Eelas = -Ebc-Ge + ^Ge-Si + -^Si-Su 

where -Esi-Si is the energy due to the interaction between the Silicon atoms. The other 
terms are analogously defined. One has 

-E-Si-Si = ^ 2^ bond ^ bond '' 2 ' 

Si— Si bonds 

where fc DO nd is either or kp depending on whether the bond is longitudinal or 
diagonal, ^bond is the amount the bond has been stretched from the equilibrium con- 
figuration. This can be written in terms of the displacement field as 

£si-Si = --uf^siusi, (27) 

where we denote by ^4si the (infinite dimensional) matrix that provides the force on all 
Silicon atoms as a function of the position of the Silicon atoms. 

The energy due to the interaction of the Germanium atoms can be written as 

-Ebc-Ge + Ece-Si = ^ ^ ^ bond (^bond ) 2 , (28) 
all Ge bonds 
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where /cbond is as in (|2f))) but here ^bond represents the amount the Germanium bonds 
have been stretched from their original equilibrium configuration (as opposed to the 
reference configuration). This can be written as 

-^Ge-Ge + ^Ge-Si = X e ^ k 

all atoms 
bonded to 
a Ge atom 

where egjk is the total elastic energy stored in all the bonds associated to the atom 
located at site (£,j, k). The factor ^ accounts for the double counting of the summation. 
We can write 

££jk = e-ijk + e \jk + e ljk 

with 





k L 
2 


{[Ul+ljk — 


U£jk + d x ] 2 + [ug^ijk — U£jk 


- d x r 


) 


+ 


k D 
2 


([ue+ijk+i 


— Uijk + d x ] 2 + [tt£_ijfc+l — 


U£jk — 


d x } 2 


+ 


k D 
2 


{[Ui+ljk-l 


— U£jk + d x ] 2 + [Ui-ljk-l — 


U£jk — 


d x } 2 


+ 


k D 
2 


{[ u e+ij+ik 


— U£jk + d x ] 2 + [u£-ij + \k — 


U£jk — 


d x } 2 


+ 


k D 
2 


{[ue+ij-ik 


— U£jk + d x ] 2 + [ue-ij-ik — 


U£jk — 


d x } 2 



where d x = a s — a g . In the above expression, each term in square brackets represents 
the contribution to the elastic energy by a pair of atoms. If no such pair exists then 
the term is not included. One can derive analogous expressions for e|- fc and e| jfc where 
d y = d x and d z = — a g . 

To compute the total elastic energy, the sum in Ea. l|28|) is computed directly, while 
the sum in Eq. (|26|) . which contains infinitely many terms, can be computed by the 
following argument. First, let us distinguish among the surface and the bulk atoms. Let 
us denote by — Fg; the force acting on Silicon due to the presence of the Germanium. 
They take the form 

where the dimension of the vector F s is equal to the dimension of the surface of the 
substrate, while represent the force acting on rest of all the infinite atoms of the bulk. 

At equilibrium, the net force acting on all Silicon atoms is zero, therefore we may 
write 

-F Si + A Si ( Us )=0. 

V u bulk / 

Using the above relations one obtains 

We remark here that F s is the same surface force computed in the previous section 
using the discrete Fourier transform. 
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5 Time stepping approximations 

The method outlined in Section [2] is unpractically slow for the following reasons 

1. For each attempted hop, a complete elastic computation has to be performed. 
Since most attempts are rejected, most of the time would be spent performing 
elastic computations that are never used. 

2. Even without elastic effects, the simple rejection-based KMC described here is 
very slow. A more effective technique would be to use rejection- free KMC, in 
which all possible events are sampled according to their probability [3J E] • How- 
ever, with the inclusion of elastic effects, the implementation of rejection-free 
KMC is not so straightforward. 

Here we shall outline various approximations of the model which lead to a much 
faster code, without significantly compromising the physical fidelity. 

As mentioned above, in order to know the rate at which an atom might hop we 
must compute the change in elastic energy of the entire system with and without that 
atom present. We make the following approximation, we assume that the change in 
the elastic energy is due to the energy in the bonds that directly connect that atom. 
Therefore we have 

&-E e i as « — ^-(4>ond) 2 = egjk (30) 

bonds to 
atom p 

where (£, j, k) is the site of the pth atom. The advantage of this approach is that 
a new equilibrium configuration has to be computed only if the move is accepted. 
Another approximation involves updating the displacement field after a given number 
J of hops. Furthermore, one can verify that the change in elastic energy for atoms 
which are lightly bonded (Nf, < 5) is very small and consequently we assume it to be 
zero. 

Finally, in order to reduce the number of rejections we separate the lightly bonded 
(iV;, < 5) and the more strongly bonded (iV& > 5) atoms in our implementation of 
kinetic Monte Carlo. This is done as follows. We take Q steps where we update the 
lightly bonded atoms and then take one step where the strongly bonded atoms are 
allowed to move. The accepted rate for the strongly bonded atoms is increased by a 
factor Q. 

The above discussion can be conveniently summarized as the following algorithm. 

ALGORITHM 

1. Choose a site at random among all M x M sites (only atoms on the surface are 
allowed to move) 

2. If there is Germanium atom present compute iV& (number of bonds) 

3. Let it hop without any elastic computation if JV& < 5 

4. If Nb > 5 ignore the atom 

5. Repeat Steps 1 to 4 M 2 times 

6. Repeat Steps 1 to 5 Q times 

7. Let the system come to mechanical equilibrium (perform an elastic solve). 
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8. Choose a site at random 

9. If there is Germanium atom present compute Nf, 

10. If N b < 5 do nothing 

11. If Nb > 5 then compute a random number, r £ U[0, 1] 

12. If r < Qexp (— 7(JVj, — 5) + A.E S ) perform a hop 

13. After J hops update the displacement field 

14. Repeat Steps 9 to 13 M 2 times. 

15. We have now advanced Q steps 

6 Numerical results 

At the present time the algorithm is to too slow to perform computations for realistic 
values for the parameter values and see effects due to elastic strain. For example, 
quantum dots observed in experiments are on the order of 20 nm. This would suggest 
that we should be computing on domains on the order of 1000 x 1000. At the present 
time the largest domain for which we can simulate in a reasonable time is 64 x 64. Since 
elastic phenomena are a bulk effect, then we must increase the spring constants to be 
unphysically large in order to observe significant elastic interaction. In our simulations 
we chose kj_, = 500 and ko = 250. These values are approximately 10 times larger 
than physical values. In addition we chose F = 10 -5 . Since the hopping rate for an 
adatom is unity then it follows that the diffusion coefficient is D = \. Therefore in 
our simulations we have D/F = 2.5 x 10 4 . From an experimental point of view this 
is small. Realistic values for D/F are typically larger than 10 5 . We have chosen high 
deposition rate to reduce the simulation time. We have taken 7 = 2. Finally we have 
used Q = 5 and J = 8. Numerical experiments revealed taking smaller values for Q 
and J did not change the answer appreciably. 

In FiguresEJEl andEl we present computations using the parameter values discussed 
above but we allow the misfit to vary. Figure shows the growth for e = 0.02. Here 
one observes layer-by-layer growth. Figure 03 presents the results when e = 0.04. One 
observes that in the initial stages of growth the morophology is similar to layer-by- 
layer growth but three dimensional islands form by nucleation type events and by the 
formation of trenches The result for the case e = 0.06 are shown in Figure E3 In 
this situation three dimensional islands form very quickly by nucleation. 

In both cases, e = 0.04 and e = 0.06 we observe Volmer- Weber growth. Sim- 
ulations were performed over a wide range of parameter values and they always re- 
vealed a sharp transition between layer-by-layer growth and Volmer- Weber growth; 
Stranski-Krastanov growth was not observed. Our results are consistent with previous 
simulations in this regard. 

7 Summary 

A numerical method for the computation of heteroepitaxial growth in the presence of 
strain using kinetic Monte Carlo has been presented. A solid-on-solid model is used 
and the elastic effects are incorporated using a linear ball and spring model. It is 
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Figure 4: Heteroepitaxial simulations with e = 0.02, all other parameter values are given in 
the text, (a) 0.5 monolayers, (b) 1.5 monolayers, (c) 2.5 monolayers, and (d) 3.5 monolayers. 
The number of gray levels is equal to the maximum height. Note: higer resolution versions 
of the figures can be found at www.math.lsa.umich.edu/~psmereka 
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Figure 5: Heteroepitaxial simulations with e = 0.04, all other parameter values are given in 
the text, (a) 0.5 monolayers, (b) 1.5 monolayers, (c) 2.5 monolayers, and (d) 3.5 monolayers. 
The number of gray levels is equal to the maximum height. 
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Figure 6: Heteroepitaxial simulations with e = 0.06, all other parameter values are given in 
the text, (a) 0.5 monolayers, (b) 1.5 monolayers, (c) 2.5 monolayers, and (d) 3.5 monolayers. 
The number of gray levels is equal to the maximum height. 
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assumed that the film is mechanical equilibrium. The strain field in the substrate is 
computed by an exact solution which is efficiently evaluated using the fast Fourier 
transform. The strain field in the growing film is computed directly. The resulting 
coupled system is solved iteratively using the conjugate gradient method. Finally we 
introduce various approximations in the implementation of the KMC to improve the 
computation speed. Numerical results show that layer-by-layer growth is unstable if 
the misfit is large enough resulting in the formation of three dimensional islands. Our 
results are in agreement with previous studies E3 EEH El El EEB] • 

Currently, we are in the process of solving the elastic equations for the deposited 
atoms using multigrid and then coupling the multigrid solver to the exact solution in 
the substrate. In addition we plan to extend the model to allow for the deposition of 
several different atomistic species. 
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